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Abstract 

A vortex-lattice method for wing aerodynamics that uses nonlinear air¬ 
foil data is presented. Two applications of this procedure are presented: 
Direct Design of a Flying Wing and Inverse Identification from wind tun¬ 
nel measurements with low-aspect ratio wings. A Newton method is em¬ 
ployed, which not only allows very fast solutions to the nonlinear equations 
but enables the calculation of static and dynamic stability and control 
derivatives without further cost. 


1 Introduction 

Fast methods to obtain aerodynamic characteristics of aircraft are and have been 
of great interest in all aircraft disciplines. The governing Navier-Stokes equa¬ 
tions result in far too complex calculations during the initial phases of aircraft 
design, and for small projects are infeasible even for post-design and evaluation. 
For the purpose of identifying the aerodynamic characteristics of aircraft, the 
tools used today are still based on the theoretical framework of the 40s. These 
are the Vortex-Step, Vortex-Lattice and Panel methods. They are all based on 
the same idea: Assuming non-rotational, incompressible flow and discretization 
this flow by combinations of fundamental solutions at certain places on the wing 
and body. The differences in these methods lies solely on the amount, and po¬ 
sition, of these solutions. This results in surprisingly accurate predictions even 
for simple discretizations used in the 40s. Today, Panel Methods are among the 
most widely used aircraft design tools. 

The drawback of these methods is simple. For Vortex-Step and Vortex- 
Lattice methods, no knowledge whatsoever is used on the shape of the airfoil, 
besides partly camber which effectively results in a change of zero-lift angle of 
attack. This results in significant errors for configurations with unconventional 
airfoils, for example high-lift configurations and low-speed small aircraft. No 
stalling characteristics are used, and no profile drag can be estimated. Even 
more detailed Euler Equation solvers cannot properly model the drag, although 
some promising developments are being made 121 HD- These approaches are still 
quite experimental though and far from being used in common applications. 

* m.ranneberg@tu-berlin.de 


1 



However, two-dimensional data of profile lift, drag and moment are available. 
Reliable simulation of 2D profiles has been feasible for years, detailed measure¬ 
ments are available and today, full Navier-Stokes solutions of the 2D problem 
are possible within hours on a common machine. The software XFLB0, if panel 
or vortex-lattice method is chosen, uses the 2D data calculated from XFOIL|T] 
to append viscous drag related to the current local lift and Reynolds-Number 
and to estimate the maximum local lift coefficient. 

Methods that try to use the 2D data in a more involved way exist since the 
40s. In [2] a method of incorporating the local lift coefficient into a lifting-line 
method is proposed. There, the nonlinear equations for the vortex strength 
along the span is solved using fixed-point iteration. This methodology is, for 
example, implemented in the software miarex , which is now part of XFLR if the 
method is chosen. The drawback are the drawbacks of the lifting-line method: 
For swept and low aspect ratio wings the results should be handled with caution. 

Another method with a similar approach was proposed in U . It uses the 2D 
data together with a generalized Vortex-Step or Weissinger method, which takes 
into account sweep and low aspect ratio. It is based on the generalization of the 
tangential-flow condition commonly applied to Vortex-Lattice Methods which 
will be discussed later. There, too, a fixed-point iteration is used for the resulting 
equations. In [5] the method has been applied to obtain simplified nonlinear 
dynamics equations for aircraft for control analysis. More recent studies, for 
example in [EJ a similar technique is used to obtain corrections to the classic 
linear Vortex-Lattice Methods. Instead of using a different boundary condition, 
an ad-hoc change in the local angle of attack is assumed. However, one could 
argue that the generalized boundary conditions of [3] lead to a similar approach. 
These conditions effectively modify the classic local downwash with another 
downwash term, resulting in a local change in angle of attack (albeit slightly 
more motivated by theory and dependent on the local vortex strength). Another 
approach been presented in [7J, which gives an excellent overview of the results 
thus far and where the change in local angle of attack is interpreted as a change 
in camber due to the flow separating from the wing earlier. There, the local 
lift and moment data is used for changing a vortex-lattice method with two 
cord-wise panels with a two-parameter decambering function and the resulting 
nonlinear equations are solved using a Newton method. Camber of an airfoil 
describes the asymmetry of an airfoil and results in non-zero lift at zero angle of 
attack. Decambering can be interpreted as a means to incorporate the increase 
in boundary layer thickness towards the trailing edge and the effective reduction 
of lift due to the reduction of (positive) camber into a three-dimensional method 
gj. Subsequent applications of this work have been published, for example in 
[5], where an approximation of the method is used for simulation and control. 
However, all methods share the same goal: Obtain a solution at every section, 
such that the global lift due to vortex strength and the resulting effective angle 
of attack is in accordance to the lift coefficient due to this angle of attack. 

Here, a method is described based on the principles derived in [|], which is 
similar to the method described in [5] . The main reason for this basis is the 
elegance of generalizing the boundary conditions beyond linear airfoils. The 
boundary conditions in [5] will be motivated by revisiting the | cord theorem 
and the relation to the generalized boundary conditions. In contrast to the 
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method employed in [7], the generalized boundary conditions do not lead to the 
generalization of two cord-wise parameter variations and thus only to a single 
cord-wise lattice discretization. 

A Newton method is employed and results in fast results with usually only 
one or two iterations during an angle sweep. The results are used for detailed 
aerodynamic analysis for static and dynamic coefficients. Derivations of these 
coefficients are possible that do not necessitate a solution to the nonlinear equa¬ 
tions, but can be obtained using the gradient of the equations. Since the gradient 
needs to be evaluated for the Newton method anyway, the dynamic derivatives 
are a convenient byproduct of the method. 

Additionally, the method is used in an inverse fashion to correct windtunnel 
results obtained with wings. 


2 The Solution Method 

The Vortex-Lattice method relies on the discretization of the solution to the 
incompressible Euler equations. Any closed path of constant vortex-strength is 
a solution to the equations, and the closed paths used here are classic Horseshoe 
paths with a tip on the quarter cord of the local wing surface, two legs along 
the cord, two legs from the trailing edge aligned with the free-stream airspeed 
and a connecting line at infinity, which results in no downwash. First, the 
basics of these methods is explained, which are the Biot-Savart Law for the 
downwash, the Kutta-Joukowski Law for the lift and the Pistolesi Theorem for 
the boundary conditions. Special attention is necessary during the formulation 
of the boundary conditions, since they are modified to enable nonlinear airfoil 
data. In the remaining section, it will always be assumed that v a is the unit 
vector of the airspeed direction. That is, |r> a | = 1. 


2.1 Downwash 


The downwash induced at some point m from a vortex path S of constant vortex 
strength T is given by the Biot-Savart law 

D m = f f ell x (1) 

47 t J 5 |(r(s) — m) | 

If the path S' is a linear path between two points ro,ri the downwash has a 
simple analytic anti-derivative, given by 


D 2’ r2 = 4 ^(n -r 0 ) 


/ r\ — m 
Vln -m\ 


r 0 -m \ (r 0 - to) x (rq - r 0 ) 

ho -m\J |(r 0 - to) x (ri - r 0 )| 2 ‘ 


( 2 ) 


All parts of the path given by a Horseshoe vortex can be calculated. Points at 
an infinite distance can be calculated, too, by using the limit as r — > oo. For 
example, the trailing leg in the airspeed direction from the trailing edge r* can 
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be written with tq = r t , n = r t + sv a and with the assumption that |i; a | = 1 as 

n r t ,oo _ y _£_/ t f r t + svg-m _ r t — m \ (r t - m) x (sv a ) 

47 r^ SVa {jr t + sv a - mj \r t -m\) \(r t - m) x {sv a )\ 2 

li m — ( W T f Tt + SVa ~ m _ n-m \ (rt - to) x v a 
s->oo 47 t 1 “ \\r t + sv a - m\ \r t -m\) \(r t - m) x v a \ 2 

T r t — m \ (r t - to) x v a 
47 T V Va \r t - m\) \(r t - m) x v a \ 2 ' 

Another special case which will be discussed during Pistolesi’s Theorem is ro = 
— sv, r\ = sv which results in 


D 


— 00,00 
m 


lim —( 2 u) T 

s—>• oo 47T 

r m x v 
2n |to x w| 2 


sv — m 
| si; — m| 


—si; — m \ (—si; — m) x (2v) 

| — sv — m\J \(—sv — m)x(2v)\ 2 


(3) 

(4) 


All of these downwash velocities are vectors. The downwash in the direction of 
interest, usually downwards ( 0 , 0 ,— 1 ) T with respect to the local surface axis, 
can be evaluated by using the scalar product with the direction of interest. 


2.2 Boundary Conditions 

With the discretized Vortex paths the possible solutions are constrained to all 
linear combinations of the individual Tj. The lift associated with the vortex 
combination is given by the Kutta-Joukowski law, which states 

L = v a x T. (5) 

The values of Tj are found by defining sufficient boundary conditions. The 
classic Weissinger method uses the tangential flow condition at the |cord. At 
this point, the three-dimensional downwash D 3 induced by all elements should 
be equal to the geometric angle of attack, resulting in an effective angle of 
attack of zero and a flow tangential to the lifting surface. However, there are 
more subtleties associated with this boundary condition. 

The theorem of Pistolesi states that at the |cord the angle of attack results 
in the correct lift associated with the vorticity. A 2D airfoil, or infinitely long 
wing, with a constant sweep 7 is considered. The lift slope of such a wing is 
given by c^ D (a) = 27r cos 2 7 a as derived from simple swept wing theory. Which, 
basically, reduces the effective airspeed due to the normal direction with respect 
to the nose by cos 7 . Since the lift depends quadratically on the airspeed lift 
is reduced by cos 2 7 . At the |cord a single vortex path with strength T runs 
along the wingspan to infinity. The lift coefficient of this wing with this vortex 
strength is given by v a x T = cosyT. The downwash at the ^ line is given by 

it = (sin 7 ,cos 7 , 0 ) T , dl = udy, /i=^-, 0,0 

D 2B = L r dl x (y« - h ) 

* 4tt J_ 00 I yu - h\ 3 

m r 

2 - 7 T cos 7 


( 6 ) 

(7) 

( 8 ) 
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Thus, the downwash results in the angle of attack associated with the vortex 
strength. 

The boundary condition of the classic Weissinger method enforces this con¬ 
dition for the three dimensional downwash. Note that the tangential-flow condi¬ 
tion uses the sweep implicitly, since the assumption is that the terms cancel each 
other out. The boundary condition of the Weissinger method can be written in 
a different way that does not use the assumption of a two-dimensional lift slope 
of 27 tcos 2 7 but instead uses a general two-dimensional lift: 

c? D (a-£)| D r + Di D r)=|u a xr| (9) 

4 4 

That is, an effective local angle of attack is used to search for a solution T where 
the effective local angle of attack results in local lift coefficients equal to the lift 
given by T. For a linear lift slope, the equation results in the tangent-flow 
conditions as already shown in [3], since the 2D downwash times the lift slope 
equals the term \v a x T| as shown above. It should be noted that the tangent- 
flow conditions are in fact a misnomer. For all other positions on the wing, for 
example the trailing edge, there is no tangential flow. It just happens to fit 
the interpretation of intuitive boundary conditions, if the |cord control point is 
chosen. 

2.3 Nonlinear Coupling 

In the previous section the boundary condition of the Weissinger method were 
reformulated without explicitly stating the 2D lift slope. These generalized 
boundary conditions © can be used with different two-dimensional lift slopes. 
And since nonlinear equations are used, it seems appropriate to use the following 
equation 

cf B (a - atan(£>l D r + Zbfr)) = v a x T, (10) 

4 4 

since the induced downwash is in fact not an angle, but the ratio of induced 
velocity and free-stream velocity which is similar only for small downwash ve¬ 
locities w.r.t. the free-stream velocity. A damped Newton method is employed 
to solve the nonlinear equations. 

2.4 Induced Drag 

The induced drag is calculated by using the common Trefftz-Plane Analysis 
far behind the wing. There, parallel to the trailing edge, the downwash is 
calculated by the same formula derived from the Biot-Savart law, but special 
care is necessary for the limit as tq —> oo. Here, the control-point m lies in the 
Trefftz-Plane and the only downwash contributions are the trailing legs. 

2.5 Washout and Angle of Attack 

In the implementation, the angle of attack is derived completely from the given 
geometry and the current airspeed vector. The washout of the airfoil sections 
is directly applied to the cord-wise legs of the vortices. This leads to non- 
symmetric legs. An example of the geometry is given in FigQ] where the icord 
line is given but the control points are not shown. 
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2.6 Geometric Angle of Attack 

For simple Horseshoe-Vortex geometries with legs and assumed airflow parallel 
to the x-axis, the geometric angle of attack is given by the angle of attack of the 
airflow and additionally the angle of attack resulting from the washout. More 
complex geometries with rotations and non-parallel legs result in less trivial 
geometric angles of attack. The surface is given by the vectors 9 = ur yz and 
the cord direction c. The local geometric angle of attack, perpendicular to this 
surface, can be calculated as follows: The goal is to find an angle ao about 
which the cord direction is rotated to result in a maximal scalar product with 
the airspeed v a . 


ao = max a (Rot ^{a)c) T v a 


( 11 ) 

( 12 ) 

(13) 


= max a (ccosa + sina0 x c + (1 — cos a)c T 99) T v a 
=>- 0 = ((—c + c T 98 ) sin a 0 + (0 x c) cos ao) T v a 



(14) 


The denominator within the atan function is only zero, if and only if the airspeed 
is completely perpendicular to the cord-wise direction. 


2.7 Sideslip 


The influence of sideslip is incorporated directly by the direction of v a and results 
in a different three-dimensional downwash due to the trailing wake as well as 
in the direction of the local lift. The geometry is unchanged, the bound vortex 
and the legs on the wing surface are still in sweep-wise direction or cord-wise 
direction, respectively. This approach of calculating the influence of sideslip has 
been originally proposed by WeissingerJTDj and is in widespread use today, for 
example in the software XFLR. It differs, however, from the approach in 0] 
where sideslip is modeled using asymmetric sweep angles. 

2.8 Stationary Coefficients 

The stationary coefficients can be deduced directly from the solution of the 
Vortex strength. With the generalized Kutta-Joukowski law, the forces due to 
lift can be calculated, and the drag is given as well. Assuming that all forces 
act on the ^cord point, the aerodynamic moments are given as well. Here, the 
airfoil moment is added to the total moment. 

2.9 Dynamic Coefficients 

The dynamic linear coefficients can be calculated without reevaluating the down- 
wash or search for a new T. Exemplary the roll-damping coefficient at some 
angle of attack with zero sideslip is derived. The roll damping coefficient is 
given w.r.t. the non-dimensional roll-rate p = Pre ° imac . If the stationary local 
velocity is vo = Voo A with v a = (cos a, 0,sina) T , the new velocity at a distance 
r from the reference point or rotation is given by v = vq + uj rea i x r. For rolling 
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only, this is v = Vo + Preai( 0, — r z , r y ). The new angle of attack and sideslip can 
be linearized with the atan2 derivative and is given by 


x' = atan2 ( v z , v x ) = atan2 ^sin a + cos a\ 


a — cos cx-r-p 
h 


/?' = atan2 ( v y , v x ) = atan2 cosa ^ 


/3 + cos a—^p 
h 


(15) 

(16) 

(17) 

(18) 


These changes in the angle of attack result in a local change of the angle of 
incidence and in a change of direction of the lift, but also in a slight change 
of the vorticity. Equation m can be approximated around the current T and 
a and solved for the changes AT due to An. With the current local effective 
angle of incidence ct e ff : the complete downwash operator D and the local lift 
slope c[(a e ff) the change in vortex strength is given by the linear equation 


/ v a x f 

V c K a e//) 


2 ) 


AT 


(19) 


With the new vortex strength and the new lift and drag directions, all forces 
and torques can be reevaluated by the Kutta-Joukowski Law and the induced 
drag definition. Finite differences can then be used to calculate the coefficient. 
That is, 

^ _ L(a, T) + L(a + Aa p , T + AT p ) 

P P 

The linear equation (1191) can be used for all dynamic derivatives, not only for the 
rotational rates but also for the aircraft velocities and linear a,/3 derivatives. 
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Figure 1: Geometry of a wing. The washout is directly included in the vortex geometry. 
The I cord line is hinted and defines the bound vortex position. 




a a 


Figure 2: Section data obtained with XFOIL. 


3 Results for a Flying Wing 

Here, calculations of one flying wing are presented. The wing is shown in Fig[TJ 
and is similar to delta-wing speed-kites. Winglets and main wing are discretized 
by 24 and 100 Vortex elements. A single airfoil along the wingspan is assumed, 
albeit with a different washout. The 2D curves of lift, drag and moment were 
calculated using XFOIL and can be seen in FigO 

3.1 Local Lift and Coefficients 

The local lift at zero sideslip at different angles of attack is shown in Figjo] 
Above 20° stall regions appear, where the local angle of attack is above the 
profile maximum. The general coefficients are shown in FigU Below and above 
the given angles of attack, no interpolation is possible and thus no results are 
available. The coefficient of maximum lift is a bit below 1.6, while the maximum 
section lift data is a bit above 1.6. Next to the stalled region of the coefficient 
curve, it is also of interest that the lift slope of the aircraft is not linear within 
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Local Lift and Drag 



Figure 3: Local lift coefficient. Winglets not shown. Stalled regions can be seen at 24° 
and 25°. 


the unstalled region of low angles of attack. This leads to different optimal 
angles of attack for best glide-ratio and best climb performance. 

3.2 Static Sideslip Stability 

The sideslip results in a change in airspeed direction and thus in a change in 
lift and drag direction. Additionally, the trailing vortex direction is changed. 
All coefficients are meant in terms of body-coordinates. That is, the sideforce 
points spanwise and the moments are meant in the aircraft coordinate system. 
The point of reference is set to be the center of gravity which is quite close to the 
mean quarter-cord line. The sideslip coefficients are the sideforce to to sideslip, 
the yawing moment due to sideslip and the rolling moment due to sideslip. They 
can be seen in Fig|5j where the coefficients were evaluated at different sideslip 
angles. However, the coefficients divided by the sideslip angle are in more or 
less the same for all sideslip angles and a linear coefficient for all sideslip angles 
can be assumed. It can be seen that the directional stability decreases with the 
angle of attack and the aircraft becomes unstable above 12°. Lateral stability 
increases with angle of attack. 


4 Inverse Identification of Airfoil Section Data 

The results so far were concerned with the estimation of aerodynamic proper¬ 
ties of 3-dimensional wings using 2-dinrensional results from airfoil calculations. 
Just as interesting is the problem of identifying the properties of the airfoil sec¬ 
tion using windtunnel measurements with 3-dimensional wings. In fact this is 
quite common, even if the wing is simply a rectangle. Instead of using end- 
plates on rectangular models or using some formulas for downwash effects, the 
method presented here can be used to estimate the airfoil characteristics that 
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Figure f: Basic aerodynamic wing coefficients. 



Figure 5: Lateral wing coefficients, evaluated at different sideslip angles. 
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lead to measurements with 3-dimensional models. The results in this chapter 
were obtained with the following procedure. Let F m = [CL m , C-Dm, (7M m ] the 
coefficients of a three-dimensional wing as measured in a wind-tunnel. The two- 
dimensional section characteristics are modeled with a linear function over the 
angle of attack with the unknown coefficients x = [cl, cd, cm] and the resulting 3- 
dimensional characteristics from the nonlinear Vortex-Lattice method is F 3 d (a;). 
Then, the solution to the following problem should result in 2-dimensional air¬ 
foil characteristics that best describe the measurements under the assumption 
that the nonlinear Vortex-Lattice method is applicable here. 

min || F 3 D (x) - F m \\ 


However, there are some problems with simply minimizing this function. The 
method calculates apparent angles of attack at every section, which are based 
on the global lift distribution. This relationship is not necessarily unique. Ad¬ 
ditionally, it can be ill-defined. Due to downwash, the maximal lift coefficients 
of the airfoil data are only occurring in the center of a wing, and the local 
angles of attack are all shifted towards the zero-lift angle of attack. For low 
aspect-ratio wings, this problem becomes more pronounced. To guide the so¬ 
lution towards sensible scenarios, a regularization penalty is introduced. Here, 
the second derivative of the airfoil data with respect to angle of attack is used 
which results in the following problem with the regularization parameter /i 

ii 112 

min \\F 3 D (x) - F m \\ + 

X 

For linear coefficient functions, this results in no penalty. The problem is solved 
by a Gauss-Newton Iteration, and the gradient of F 3 d with respect to x is 
computed using finite differences. The method has been applied to the identi¬ 
fication of the Selig S1223 [12] profile from a rectangular wing with an aspect 
ratio of 2.7. The airfoil exhibits highly nonlinear behavior and is prominent in 
low Reynolds-Number applications. Good measurements for this airfoil have 
been obtained at U1UC [~HZ1 . with which the measurements and reconstructions 
will be compared. 

The measurements were taken at a small windtunncl at the BIONIK Insti¬ 
tute at the TU Berlin. In FigJT] the measurements taken in Berlin, together 
with the artificial measurements obtained using the presented method with the 
reconstructed polars as well as the polars from the UIUC. In FigjHlthe 2D Lift 
coefficients reconstructed with the inverse method from the wing measurements 
and the artifical data can be seen, as well as the original S1223 measurements 
at UIUC. 

At Lift coefficients between 1 and 1.5, the reconstruction agrees well with 
the measurements of the UIUCQ21. At the Reynolds Number the maximum 
Lift coefficient should exceed 2.0, which is not reconstructed here. Instead, the 
section lift seems to be stalling earlier. Additionally, the characteristic drop in 
lift below 2° is not reconstructed but only hinted. 

Using the artificial data, the 2D profiles are reconstructed very well un¬ 
til maximum lift, only the characteristic drop at around -2° has been slightly 
smoothed. Hence, the inverse calculation works rather well. However, using the 
measurements taken in Berlin, the airfoil data is not reconstructed well. Either 
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the wing has not been constructed close enough to the original S1223 specifica¬ 
tions or the nonlinear Weissinger method is not applicable to these low aspect 
ratios. 
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2.5 



Figure 6: Comparison of Section Lift from the SI223 airfoil reconstructed from mea¬ 
surements at the windtunnel at the BIONIK Institute at the TU Berlin with Measure¬ 
ment Data from the UIUC U2t . AR 2.7, RE 200,000. 



Figure 7: Comparison of Wing Lift with the SI223 airfoil at the windtunnel at the 
BIONIK Institute at the TU Berlin with measurements and forward calculation using 
Measurements from the UIU C\12t . AR 2.7, RE 200,000. 
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